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Abstract 

Many standard approaches to geometric model fitting are based on pre- 
matched image features. Typically, such pre-matching uses only feature ap- 
pearances (e.g. SIFT) and a large number of non-unique features must be 
discarded in order to control the false positive rate. In contrast, we solve fea- 
ture matching and multi-model fitting problems in a joint optimization frame- 
work. This paper proposes several fit-&-match energy formulations based on 
a generalization of the assignment problem. To find near optimal solutions, 
we developed an efficient solver based on min-cost-max-flow algorithm. Our 
approach significantly increases the number of detected matches. In practice, 
energy-based joint fitting & matching allows to increase the distance between 
view-points previously restricted by robustness of local SIFT-matching and to 
improve the model fitting accuracy when compared to the state-of-the-art. 

1 Introduction 

Many existing methods for model fitting and 3D structure estimation use prematched 
image features as an input (bundle adjustment pQ, homography fitting [21 [3], rigid 
motion estimation [U EJ E]). Vice versa, many matching methods (sparse/dense 
stereo) often use some preestimated structural constraints, e.g. epipolar geometry, to 
identify correct matches/inliers. This paper introduces a novel framework for simul- 
taneous estimation of high-level structures (multi-model fitting) and low-level corre- 
spondences (feature matching). We discuss several regularization-based formulations 
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of the proposed fit & match (FM) problem. These formulations use a generalization 
of the assignment problem and we use efficient specialized min-cost-max-flow solver 
that has been overlooked in the computer vision community. This paper primarily 
focuses on jointly solving multi-homography fitting and sparse feature matching as 
a simple show case for the FM paradigm. Other applications would be rigid motions 
estimation, camera pose estimation [7J, etc. 

Related Work: In case of reliable matching, RANSAC is a well-known robust 
method for single model fitting. The main idea is to generate a number of model 
proposals by randomly sampling the matches and then select one model with the 
largest set of inliers (a.k.a. consensus set) with respect to some fixed threshold. In 
case of unreliable matching, e.g. repetitive texture or wide view-point, RANSAC or 
any technique that relies on pre-computed matching would fail. 

Guided-MLESAC [8j and PROSAC are RANSAC generalizations that try to 
overcome unreliable matches while generating model hypotheses. Their main idea 
is to ensure that matches with high matching scores are more likely to get sampled, 
thus "guiding" the sampling process while generating model hypotheses. One could 
argue that these techniques would still fail since false matches could also have high 
matching scores, e.g. scenes with repetitive texture. SCRAMSAC [TU] is a form of 
spatial guided sampling that uses a spatial consistency filter to restrict the sampling 
domain to matches with similar local geometric consistency. This method is sensitive 
to the ratio of occluded/unoccluded features, as in that case the assumption that 
correct matches form a dense cluster is no longer valid. The main drawback of these 
RANSAC generalizations is that they focus on generating a reliable model hypotheses 
by using pre-matched features (fixed matching). That drawback could be avoided 
by jointly solving the matching and fitting problems. 

An attempt to formulate an objective function for fitting- & -matching naturally 
leads to a version of the assignment problem. The majority of prior work could 
be divided into two major groups: matching techniques using quadratic assignment 
problems and FM techniques using linear assignment subproblems. 

Quadratic assignment problem (QAP) normally appears in the context of non- 
parametric matching. For example, the methods in (TTJ [121 03] estimate non-rigid 
motion correspondences as a sparse vector field. They rely on a quadratic term in 
the objective function to encourage geometric regularity between identified matched 
pairs. Such QAP formulations often appear in shape matching and object recogni- 
tion. QAP is NP-hard and these methods use different techniques to approximate it. 
For example, [H] approximates QAP by iteratively minimizing its first-order Taylor 
expansion, which reduces to a linear assignment problem (LAP). 
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If correspondences are constrained by some parametric model(s), matching often 
simplifies to LAP when model parameters are fixed. In this case, the geometric 
regularity is enforced by a model fidelity term (linear w.r.t. matching variables) and 
pair-wise consistencies [TT| [T2l [13] are no longer needed. Typically for FM problems, 
LAP-based feature matching and model parameter fitting are preformed in a block 
coordinate descent (BCD) fashion. For example, SoftPOSIT [7] matches 2D image 
features to 3D object points and estimate camera pose in such iterative fashion. 
Building on the ideas in SoftPOSIT Serradell et al. [15] fit a single homography using 
geometric and appearance priors with unknown correspondences. It should be noted 
that although SoftPOSIT solves for the matching by using softassign [14J — which will 
result in an approximate matching. SoftPOSIT justified using softassign by claiming 
that the objective function is non-linear. That claim is not exactly accurate as 
SoftPOSIT uses BCD to minimize that non-linear objective function. Thus, when 
solving for the matching all the other parameters are fixed and the objective function 
becomes linear in terms of the matching variables. For linear objective function the 
optimal matching could be found by using a Linear Programming solveiQ 

Our work develops a generalization of LAP for solving FM problem when match- 
ing is constrained by an unknown number of geometric models. In contrast to [3 US], 
we do not assume that matches/correspondences are constrained by a single paramet- 
ric model. Note that in order to solve FM problem for multi-models, a regularization 
term is required to avoid over fitting. Unlike [T5l [7] , our energy formulation includes 
label cost regularization [TB] and spatial regularization as in some QAP approaches 
to matching 

Contribution: In this work we propose two FM energy functionals ^ and ^ for 
jointly solving matching and multi-model fitting. Energy ^ consists of two terms: 
unary potentials for matching similar features and assigning matched features to 
their best fitting geometric models, and a label cost term to discourage overfitting 
by penalizing the number of labels/models assigned to matches. Energy ^ consists 
of unary potential and label cost terms, as in energy @, and a pairwise potential 
term for encouraging nearby matches to be assigned to the same label/model. 

The key subproblem when minimizing ^ or ^ in BCD fashion is generalized 
assignment problem (GAP), which is our novel generalization of LAP to mutli-model 
case. GAP jointly formulates feature-to-feature matching and match-to-model as- 
signment. We propose a fast globally optimal approach for GAP by generalizing 
min- cost-max-flow techniques for single model LAP matching [TT] . 



1 SoftPOST use of softassign to find the matching could be attributed to performance gains, as 
LP solvers will be considerably slow for large problems. 
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Figure [8] compares the results of a standard energy-based multi-model fitting 
algorithm [16] (EF) and our proposed energy-based multi-model fitting- & -matching 
algorithm (EFM). EF used the standard pre-matching technique in [18] that rejected 
a relatively large number of true matches. EFM found better models' estimates 
because it nearly doubled the number of identified matches. 

2 Our Approach 

Standard techniques for sparse feature matching p2] independently decide each 
match relying on the discriminative power of the used feature descriptor. These tech- 
niques are prone to ignoring a large number of non-distinct image features that could 
have been valid matches. Our unified framework simultaneously estimate high-level 
structures (multi-model fitting) and low-level correspondences (features matching). 
Unlike standard techniques, our approach is less vulnerable to the descriptor's dis- 
criminative power. We discuss regularization-based formulation of the proposed fit 
& match problem. While there are many different applications for a general FM 
paradigm, this work primarily focuses on jointly solving geometric multi-model fit- 
ting (homographies) and sparse feature matching. 

We will use the following notations in defining our energy: 



J-~i - set of all observed features in the left image. 

T r - set of all observed features in the right image. 

C - a set of randomly sampled homographies (labels). 

f p - label assigned to feature p such that f p G C 

f - a labeling of all features in the left image, / = {f p \p G F{\ 

dfi - parameters of homography h from left image to right image. 

9 - set of all models' (homographies) parameters. 

Si - Subset of features in the left image supporting one geometric model 

(plane, homography), see Figs Q(a-b). 
S r - Subset of features in the right image supporting one geometric model 

(plane, homography), see Figs (ij)(a-b). 

binary variable which is 1 ifp and q are matched (assigned) 

to each other and otherwise. 
•M := {x pq | (p, q) GJfX J>}. 



Q(p, q) - appearance penalty for features p G T\ and q G T r 
based on similarity of their descriptors. 
Af - edges of near-neighbor graph, e.g. delaunay triangulation, 
for left image features. 
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2.1 Energy 

We will define the overall matching score between two features p G J 7 / and q e J> as 
a function of geometric transformation 6^ 

D pq (6 h ) = \\0 h -p-q\\+Q(p,q) (1) 

combining the geometric error and the appearance penalty where 1 1 1 1 denotes geo- 
metric transfer error. A similar matching score was used in computing the ground 

truth matching in [THl 120] • We can also use a symmetric matching score 

D m {9 h ) = \\6 h -p-q\\ + \\0? • q -p\\ + Q(p,q). (2) 

We are only interested in symmetric appearance penalty Q(p,q), e.g. the angle (or 
some metric distance) between the features' descriptors of p and q. From here on 
D pq refers to the symmetric matching score. 

In this work, Q(p,q) = if the angle between the two features' descriptors is 
less than 7r/4 and oo otherwise. The aforementioned non- continuous appearance 
penalty is less sensitive to the descriptor's discriminative power in comparison to the 
continuous one. 

To simplify our formulation we will introduce our energies under the assumption 
that there are no occlusions 

E l (f,9,M) = ^W^ M + ^M/) (3) 

p&Ti hec 

s-t. E ff6JV ^ = i Vpejr, K (4) 

Xpq e {0,1} Vpe T h VgG T t ) 

where 5h(f) = [3p <E Fi : f p = h] and [•] are Iverson brackets, and 



under contraints Q. We will show how to handel outliers/occlusions in Section 3.1 



Ei is more powerful than E\ because the spatial regularizer eliminates the artifacts 
that results from using only one regularizer in E\. The reader is refered to [16] for a 
more detailed discussion comparing E\ and E 2 for fixed matching in the context of 
multi-model fitting. 
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2.2 Optimization 

In this section, we describe an efficient approach, EFMi, to minimize E\ in a block 
coordinate decent (BSD) fashion, and a second approach, EFM 2 , to minimize E 2 . 



Energy-based Fitting & Matching for E\ (EFMi) 

Initialization: Find an initial AI using standard matching technqiues 
repeat 



Given At, solve (m using PEaRL |16j for / and 6 



Given 6, solve (|7j)-([8j) using LS-GAP, see Sec. 
until E\ converges 



3.2 



for M and / 



EFMi finds an initial matching using standard matching techniques and then it 
iteratively minimizes E\ by alternating between solving for / and 9 while fixing Ai, 
and solving for / and Ai while fixing 9. Although EFMi is guaranteed to converge 
since E\ is bounded below, i.e. E\ > (5, it is not trivial to derive a theoretical bound 
on the covergance rate and approximation ratio of EFMi. However, in Section [5j we 
emprically show that EFMi converges in a few iterations to a near optimal solution. 

On the one hand, E\ for fixed Ai reduces to 

peTi hec 

where D p {6h) = D pq (8h) V/i G £ provided that q is assigned to p by Ai, i.e. x pq = 1. 
Furthermore, energy (JgJ) could be effeciently solved for /, 9 using PEaRL [IB] . 

On the other hand, E\ for fixed 9 could be re-parameterized and written in the 
following form 

E(M f )= E E Y,D m {9 h )x pqh + PY, h zMM f ) (7) 
hec peFi q^T r 

hec peTi 

s-t. E Expqh = l VpGJi } (8) 

h&C q&Tr 

x pqh G {0, 1} Vh <E C, p<E Ti, q G T r J 

where binary variable x^q^ is 1 if p and g are matched to each other and assigned 
to model h, and otherwise. Matching Aif is defined as {xp^h | (p, q, h) G Ji x 
J> x £} encapsulating information of both feature-to-feature and match-to-model 
assignments, and ^ is now defined as 5h(Aif) = [Bp G Jj, g G J> : = 1]. We 
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will refer to the special unregularized case of optimization problem (j7])-(|8]) where 
(3 = as the generalized assignment problem (GAP) [^J 

This is a weighted matching problem over a fixed set of multiple models that 
match features and assigns each match to a model. GAP is an integral linear program, 
see Appendix [B] for proof, and therefore any Linear Programming toolbox could be 
used to find its optimal solution by solving its relaxed LP — but will be considerably 



slow due to the size of the problem at hand, as shown in Section 3.2 A fast approach 



to solve GAP is described in Section 3^ The optimal solution for GAP may overfit 
models to data since the number of models is not regularized when j3 = 0. For 
(3 > optimization problem ([7])-((8]) could be solved using LS-GAP, introduced in 



Section 3.2 , which utilizes a fast GAP solver in a combinatorial local search algorithm. 



This local search over different subsets of L selects a solution reducing energy (J7|). 

Energy-based Fitting & Matching for E 2 (EFM 2 ) 

Initialization: Find an initial A4,f,8 using EFMi 
repeat 

Given Ai, solve S using PEaRL for / and 6 
Given f,6, solve @-((4j) using LC-GAP, see Sec. 
until E2 converges 



3.3 



for M 



EFM2 uses EFMi result as an initial solution and then iteratively minimizes E2 
by alternating steps solving for / and 9 while fixing M., and solving for M. while 
fixing 6 and /. Energy ^ for fixed M. reduces to 

E(f, e) = e a>(%) + x E l/pt^+^E w) ( 9 ) 

and is solved using PEaRL. For fixed / and 9 energy ^ reduces to 

E(M) = E I>m(*/>m ( 10 ) 

under constraints Q. Energy (10) is solved using label constrained GAP (LC-GAP), 



Section 3^ LC-GAP is a varirant of the fast GAP solver that can change feature- 
to- feature matching without affecting the current labeling. It should be noted that, 
based on our experience, EFM2 slightly modifies the initial solution by rejecting or 
correctly matching less than a handful of false positives. Therefore, in practice we 
suggest to run only one iteration of EFM 2 to reject the false positives incorrectly 
matched due to lack of spatial coherency. 

2 Our definition of GAP is different from the formal definition in optimization literature. 
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3 Algorithms 



We introduce a fast GAP solver, in Section 3.1 to minimize energy 0-(§ for = 0. 
Section 3.2 covers LS-GAP algorithm which is used in EFM4 to minimize energy (j7j)- 
for (3 > . Finally, in Section 3.3, we describe LC-GAP which is used in EFM 2 



to minimize energy (10). 



3.1 A fast GAP solver 

As stated earlier, GAP is an integral linear program and therefore any LP solver 
could be used to find its optimal solution. However, LP solvers will be slow con- 
sidering large number of decision variables for GAP instances in our application. 
In this section, we show that our GAP problem (equations ([7])- ^ for (3 = 0) can 
be equivalently reformulated as a standard min-cost-max-flow (MCMF) problem for 
which very efficient specialized solvers are known. In the case of a single model, it is 
well-known that an assignment problem (AP) can be solved by finding MCMF on a 
particular graph [T7] . One of our contributions is a more general graph construction 
for formulating multi-model assignment problem, i.e. GAP, as an instance of MCMF. 

We will first state the formal definition of MCMF before introducing our formu- 
lation. Let Q = (V, £) denote a graph with vertices V such that s, t G V, and edges 
S where each edge (v, w) G £ has a capacity u{y , w) and cost c{v , w). The cost of on 
an arbitrary valid flow function F is defined as 

cost(F) = c(v, w) ■ F(v, w). 

A flow that maximizes the amount of flow passed from s to t while minimizing the 
cost is the MCMF. 

To formulate GAP as MCMF problem we build graph Q* = (V, £), see Fig. [TJa), 
with the set of nodes 

V ={s, t} U {n p I p e Ti} U {n q I q G T r } U 

{n p h I V G JFi, h G £} U {n qh \ q G J>, h G £}, 

the set of edges 

£= {(s,n p ) 
{(n p ,n p h) 
{(n p h,n qh ) 
{(n qh ,n q ) 

{( n q,t) 



\pefi\u 

I p G J 7 /, h G C} U 

I p G Ti, q G J>, /i G £} U 

I q G J>, /i G £} U 

I 9 e 



Tech.Rep.#755, CS dept., UWO, London, Canada, ISBN: 978-0-7714-2980-4 



and the following edge capacity u and edge cost c functions 

u{v,w) =1 for (v,w) G £ 

D pq (0 h ) for (v, w) G {{n ph , n qh ) \ p G T\, q G T r , h G £} 



c(v, w) 







otherwise. 



Lemma 3.1 The optimal solution for a feasible GAP (eq. ([7])-((8| for (3 = 0) 

Mf = {x pqh = F*(n ph , n qh ) | p G T h q G T r , h G £} 
w/iere F* : £ -> {0, 1} is the MCMF over graph Q* . 



is 



Using Lemma |3.1[ see proof Appendix |A| a GAP solution Ai / could be found by 
using an effcient MCMF algorithm |17j . 




(a) Q* of a generic GAP instance 



(b) Example of LC-MCMF Q) 



Figure 1: Figure (a) shows the generalized graph construction Q* of a generic GAP 
instance — with unit capacity edges and edge cost function c(v,w) = D pq {9h) for all 
(v,w) G {{n p h,n q h) \ Vp G Ti,q G F r ,h G C} and otherwise. This construction does 
not assume that = |J>|. Figure (b) shows G*j of a GAP with = |J>| = 3 and 
\C\ = 3 under the labeling constraint / = [1 3 2]. 
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Occlusions/Outliers: GAP might be unfeasible due to occlusions and in that case 
1.7-} | |J>|. We add — |J>|| dummy features, with a fixed matching cost T, to 
the smaller set of features to ensure GAP feasibility. This is equivalent to changing 
a rectangular assignment problem to a square one. Also, to make our approach is 
robust to outliers we introduce an outlier model <fi such that D pq (<j)) = T for any 
T\ and qE T r . The use of an outlier model with a uniformly distributed cost is a 
common technique in Computer Vision [TB| |2"T] . 



3.2 Local Search-GAP (LS-GAP) 

Now we introduce a local search algorithm that solves regularized GAP ([7])-(|8| with 



(3 > using GAP algorithm in Section 3.1 as a subprocedure. Assume that C is 



the current set of possible models^] Let C c be an arbitrary subset of C and M.f(C. 
denote the GAP solution when the label space is restricted to C c . Note that GAP 
ignores the label cost term in ([7| but we could easily evaulate energy (J7]) for A4f(C c ). 
The proposed LS-GAP algorthim greedily searches over different subsets L c for one 
such that Aif(£ c ) has the lowest value of energy ([7]). Our motivation to search 
for minima of energy ([7])-(j8]) only among GAP solutions comes from an obvious 
observation that a global minima of ([7])-(|8]) must also solve the GAP if the label 
space is restricted to a right subset of C. 

We define sets of all possible add, delete and swap combinatorial search moves as 

Af a (£ c ) = U heC \ Cc {C c Uh} 

M d {C c ) = u heCc {£c\h} 

Af s (£ c ) = U heCc {C c U£\h}. 
eec\c c 

These are three different local neighbourhoods around £ c . We also define a larger 
neighbourhood Af* around C c which is the union of the above 

N\C C ) = N a (C c )uN d (C c )uN s (C c ). 

LS-GAP uses a combination of add, delete and swap moves, as in [22J, to greedily 
find a set of labels near current set C t that is better w.r.t. energy (|7j). 



3 In practice, C is restricted to be the set of models that are assigned to at least one matched 
pair of features in energy ^ solution. 
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LS-GAP 

C t <- <f> 

while 3 C c G N t 

if energy of Mf{C c ) < energy of M f (C t ) 
C t <- C c 
Nt^N*{C t ) 
else 

Mt^M t \ C c 
return the GAP solution Aif(Ct) 



In LS-GAP, we initially set Ct to <p but it could be any arbitrary subset of C. 
Then J\f*(£ t ) is searched for a move with a GAP solution Aif(C c ) of a lower energy 
(J7|) than the current one, i.e Ai f(C t ). Such a move is accepted if it exists and 
C t is updated accordingly otherwise LS-GAP terminiates. LS-GAP will definitely 
converge since energy ^ is lower bounded and C is finite. 

The subprocedure used in solving GAP is the bottle neck of LS-GAP. Figure [2] 
compares two LS-GAP implementations; the first uses Mosek [23J to solve the relaxed 
LP of the integral linear program and the other solves GAP using the fast GAP solver 
introduced in Section 3^ Figure [2] shows that the MCMF based implementation 



(fast GAP solver) outperformed the LP one in the cases of fixed \C\ and increasing 
number of features, and fixed number of features and increasing |£|. 



3.3 Label Constrained-GAP (LC-GAP) 

LC-GAP solves a GAP instance with fixed labeling /, i.e. each left feature p must be 
assigned to a predefined model f p . LC-GAP uses a slightly different graph construc- 
tion than Q* that enforces the required labeling constraints. The graph construction 
corresponding to a GAP instance under labeling / contraints is Qj = (V, £) where 

V={s,t}U{n p \Vpe J}} 

U {n q \\/q E F r } U K/jVp G 7i} 
U {n qh \Mq eT r ,he C} 

and £, capacity function u and cost function c are as defined as in Q* provided that 
both edge nodes exist in V of Q* f , see example in Fig. [lib). 
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Label set size 



(a) LS-GAP performance for increasing (b) LS-GAP performance for increasing \C\ 
number of possible matches 

Figure 2: Figures (a-b) show the average performance, over 50 runs of randomly generated 
GAP instances, of the LP and MCMF based implemeantion of LS-GAP. Figure (a) shows 
average performance for \C\ = 5 and increasing number of possible features x |J>| where 
|J>| = \J~i\- Figure (b) shows average performance for = |J>| = 10 3 and increasing \C\. 



4 Ground truth 

The ground truth is computed by first manually identifying and segmenting regions 
corresponding to separate models (planes/homographies), see Fig. [3} Then, we com- 
pute an optimal matching of extracted features inside each identified pair of corre- 
sponding regions with respect to the geometric fitting error and appearance. This 
method is similar to the one used in [20J IT5] . 

Below we describe our technique for computing the "ground truth matching" for 
each model with manually identified spatial support, as illustrated in Fig. [3j We 
compute sets of SIFT features Si and S r inside each pair of manually identified 
corresponding regions, see Fig. |4j It is possible to independently fit one homography 
to each pair of corresponding sets {Si, S r }. For simplicity, we first assume that there 
are no occlusions, i.e if a feature appears in the left image then its corresponding 
feature appears in the right image and vise versa. Thus, the number of left image 
features equals the number of right image features. We will show how to deal with 
occlusions later. 

The SIFT features of two corresponding sets, namely Si and S r see Fig. |4j are 
matched using the criteria described in [TS]. Then we use RANSAC [21] to find a 
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homography 6h that maximizes the number of inliers between the features in two 
corresponding regions. Using RANSAC in this case is not problematic since features 
in Si and S r support only one homography/model. This homography is only used as 
an initial guess in finding the ground truth model. 




Figure 3: In this example we identified only two planes. The manually identified corre- 
sponding support regions for these two models are shown in blue and red. 





(a) Example of corresponding sets Si and S r supporting the blue model in Fig. 





3 



(b) Example of corresponding sets Si and S r supporting the red model in Fig. jij 

Figure 4: Two examples of the corresponding sets of features Si and S r supporting the 
blue (a) and red (b) models. 



Given a homography 6h, the problem of finding an optimal one-to-one matching 
that minimizes the total sum of matching scores between the left and right features 
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in two corresponding regions could be formulated as an assignment problem 

AP : argmin ]T ^ D Pq (8h) x pq (11) 

M p£Si q£S r 

s.t. ^ x pq = 1 Vg G S r 

pes, 

X pq = 1 Wp G Si 

q£S r 

Xpq G {0, 1} Vp G 5,, Vg G 5 P . 

The two linear constraints in AP enforce one-to-one correspondence between the 
features in S r and Si, see Fig. |5|a). 

For any fixed matching M. the appearance term E pgS; J2 q( zs r Q{Pi l) x pq * n AP's 
objective function becomes constant. After finding an optimal At for AP, we could 
further decrease the objective value by re-estimating homography 9h minimizing the 
geometric error, e.g. see first term in Q, over all the currently matched features. 
We can continue to iteratively re-estimate matching Ai and homography 9h until the 
objective value of AP could not be reduced anymore. 

The described optimization procedure maybe sensitive to the initial homography 
found by RANSAC. In an effort to reduce such sensitivity we repeat the whole 
procedure several times, and report as ground truth matching Ai and model Qh that 
have the lowest value of the objective function of AP. 

Now we can discuss possible occlusions that we ignored so far. The presence 
of occlusions or outliers introduces two problems. First, the number of features in 
corresponding regions are no longer guaranteed to be the same. Such an imbalance 
between the features has to be addressed in order to enforce one-to-one correspon- 
dence. Second, we can no longer assume that there exists one homography that fits 
all features. 

To balance out any possible difference between the sizes of sets S\ and S r we use 
dummy features with a constant matching cost penalty, as detailed below. Without 
loss of generality we can assume that the number of extracted features in the left 
region is less than or equal to those in the right region^} In this case, there are at 
least \S r \ — \Si\ occluded features in the right image. As illustrated in Fig. (5^b), 
an imbalance between the number of features renders the one-to-one correspondence 
constraints in AP infeasible. One way to overcome this problem is to add \S r \ — \Si\ 
dummy features to set Si, see Fig. ^c). We define a fixed matching cost penalty 
Dd q = T for assigning any dummy feature d to any q in S r . It is also possible to use 



We could always swap the regions to satisfy that assumption as long as the used geometric 
error and appearance measures are symmetric. 
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(c) Balancing with dummy features. (d) GAP with occlusion model <F 
Figure 5: Figure (a) shows the straight forward case, no occlusions, for enforcing the one- 
to-one correspondences. Notice that in this case the number of features in both images is 
the same and therefore the one-to-one correspondences constraints are balanced. Figure 

(b) shows a case with unbalanced one-to-one correspondence constraints, i.e there are no 
enough features in Si to enforce the one-to-one correspondence constraints. And, figure 

(c) shows how to balanced the one-to-one correspondence constraints by introducing the 
dummy feature d. Figure (d) shows how to account for occlusions, using occlusion model 
4>, in case of balanced one-to-one correspondence constraints. 
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only one dummy feature d in Si but for that specific feature constraint Y^ q es r Xd i = ^ 
would have to be replaced by ^2 q&Sr Xdq = \S r \ — \ S{\{d} \ . We adapt the first approach 
with multiple dummy features only to simplify our notation and avoid the special 
handling of feature d. The use of dummy feature/entity is a common technique for 
balancing out unbalanced assignment problems in operations research |25j . 

Even under the assumption that \Si\ = \S r \ occlusions are possible and we can not 
assume that there exists a homography that fits all features. In order to make our 
approach robust to occlusions/outliers we use generalized the assignment problem, 
GAP, to allow each feature to choose between two models: a homography 6h and an 
occlusion model such that D pg ((p) = T for any p G Si and q G S r . The use of an 
occlusion (or outlier) model with a uniformly distributed matching cost is a fairly 
common technique in Computer Vision [TBI E3]> see Fig- P}A). 

5 Evaluation 

In this section, we compare the matching quality of the EFM framework vs. stan- 
dard SIFT matching [18] . Then we discuss some of the EFM framework properties, 
e.g. convergence rate and the effect of the initial set of proposals size on the match- 
ing quality. Finally, we compare the quality of the estimated models by the EFM 
framework to the models estimated by an EF algorithm PEaRL[16J. 

Our matching evaluation criterion is based on Receiver Operating Characteristics 
(ROC) of the True Positive Rate (TPR) vs. the False Positive Rate (FPR). The ROC 
attributes for matching M. and ground truth (GT) matching M.gt are defined as 
follows: 

Positive (P) number of matches identified by M.gt 

Negative (N) number of potential matches that were rejected by Aior, 

i.e. N = | J} | x |J" r | - P 
True Positive (TP) number of matches identified by M. and M.qt (intersection) 
False Positive (FP) number of matches identified by M. but were rejected by M.qt 
True Positive Rate (TPR) ^ 
False Positive Rate (FPR) 

Figure |6](a) shows the ROC curve of standard SIFT matching achieved by varying 
the second best ratio (SBR) threshold where SBR is the ratio of distance between a 
left feature descriptor and the closest right features descriptor to the distance of the 
second closest. EFM is non-deterministic and the energy of convergence, a.k.a final 
energy, depends on the size of initial set of proposals \C\. Therefore, for EFM we 



Tech.Rep.#755, CS dept., UWO, London, Canada, ISBN: 978-0-7714-2980-4 



17 



show a scatter plot that relates the ROC attributes to the final energy (color coded) 
by varying \C\. As can be see, EFM outperforms standard SIFT matching and the 
lower the final energy the better the matching quality. Furthermore, Fig. |6](b) shows 
multiple histograms relating the final energy frequencies, of 50 runs, to \C\ (color 
coded). As can be seen, the bigger \£\ is the more likely the final energy is going to 
be small and the more likley that EFM behaviour becomes more deterministic over 
different runs. 

Figure [7] shows the effect of EFM iterations on the energy with respect to time 
for different \C\. For each \C\ the experiment is repeated 50 times. On the average 
each iteration took 1.5 mins, and most of the energy was reduced in the first three 
iterations. EFM converged on the average after 5 iterations. The plots in Fig. [6] 
and [7] are shown for the Metron College example, in Fig. [8j to illustrate the general 
characteristics/behavior of our method. It will be meaningless to average these plots 
over multiple examples since they would not share the same energy scale, i.e. a low 
energy for one example could be high for another one. 
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(a) EFM vs SIFT matching quality (b) Effect of \C\ on the final energy 

Figure 6: Figure (a) shows the ROC curve of the standard SIFT matches by varying the 
SBR threshold, and the scatter plot represents EFM results for different sizes of initial 
set or proposals. The scatter plot is color coded to show relation between the achived 
final energy and the quality of the matching, the lower the energy (blue) the better the 
matching. Figure (b) shows mut- liple histograms of the final energies for different sizes of 
initial set of proposals — blue indicates a large initial set of proposals while red indicates a 
small set. The larger the size of the initial set of proposals the more likely that EFM will 
converge to a low energy. 



Tech.Rep.#755, CS dept., UWO, London, Canada, ISBN: 978-0-7714-2980-4 



18 





1 .5392 




1.5382 




1.5372 




1.3611 




1.3601 


to 






1.3591 


C 




w 


1.189 




1.188 




1.187 




1.186 




1.185 



| L | =1000 
| L | =500 
| L | =100 




4 6 8 10 

Time in minutes 



16 



Figure 7: EFM energy over time in minutes. EFM converges on the average in 5 iterations, 
and an iteration on the average takes 1.5 minutes. 



Figures [8|a) and (b) show left features of EF [TH] inliers and left features of 
matches identified by EFM, respectively. Also, outliers or unmatched features are 
shown as x. EFM on the average found double the number of matches compared to 
using EF and SIFT standard matching. Figures [8]^c) and (d), and (e) and (f) are 
the zoom in for Segment 1 and Segment 2 in (a) and (b), respectively. Figures |8](g) 
and (h) show the matching, over a small region, between the left and right images 
results of EF and EFM, respectively. 

Figure [9] shows more results comparing EFM vs. EF and SIFT standard match- 
ing. In general EFM was able to find more matches than EF but EFM outperformed 
EF in two particular examples; the graphite example, shown in second row, in which 
large viewpoint between left and right images resulted in SIFT standard matching 
producing only 76 potential matches, and redbrick house example, shown in third 
row, in which repetitive texture of the bricks reduced the discriminate power of SIFT 
descriptor. 

In order to evaluate the quality of the estimated model 8^, we will use the following 
geometric error ratio GQ{6 h ) 

^. Q . STE{B h ,f GT ,M GT) 
STE(6 GT , f GT , Mgt) 

where for is the ground truth labeling and STE(8h, f,M) is the Symmetric Transfer 
Error of Oh, i.e. geometric error, computed for labeling / and matching Ai — the close 
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GQ{6 h ) is to 1 the better the model estimate. 

Table [T] shows the effect of increasing viewpoint angle, between the left and right 
images, on the quality of model estimates. As the viewpoint increases the number of 
matched points by EF sharply decreases while for EFM the decrease is not as steep. 
In addition, EF becomes more sensitive to the used SBR threshold for increasing 
viewpoint, see variance for large viewpoint. Furthermore, EFM archives near optimal 
matching, see TRP and FPR. 

The used fitting threshold T affects the ground truth, EFM and EF results, as 
it is a parameter for these methods. Table [2] shows the effect of increasing T on EF 
and EFM. For the case of T < 1, T is underestimated and running the ground truth 
multiple times will result in similar final energies but slightly different matching. The 
more we decrease T the more different the matchings will be. For T < 2 and T < 3 
the ground truth result become more deterministic over multiple runs. Finally, when 
computing the ground truth for all the examples shown above we manually handed 
tuned T to find the smallest T that gives a stable ground truth over multiple runs. 

6 Conclusions 

We introduced two energy functionals that use different regularizers for the fit-&- 
match problem. We also introduced optimization frameworks for these functionals. 
Our experimental results show that our energy-based fit- & -match framework finds 
a near optimal solution for the feature-to-feature matching and better model esti- 
mates in contrast to state-of-the-art energy-based fitting frameworks, e.g. PEaRL. 
In addition, we showed that for a given set of models it is possible to efficiently find 
the optimal feature-to-feature matching and match-to- label assignment. Our frame- 
work could be used to fit-&-match more complex models, e.g. fundamental matrices, 
without affecting the framework's complexity. It could also be used to fit-&-match 
a mixture of different models, e.g. homographies and affine transformations, and pe- 
nalize each model/label based on its complexity. Finally, we plan on applying our 
framework to estimate camera pose. 
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(a) EF left image result 




(c) Enlarged Segment 1 in (a) 




(e) Enlarged Segment 2 in (a) 




(g) Part of the EF matching, between left 
and right images, i.e. inliers of the SIFT 
standard matching. 




(d) Enlarged Segment 1 in (b) 




(f) Enlarged Segment 2 in (b) 




(h) Part of the EFM matching, between 
left and right images. 



Figure 8: Metron College from VGG Oxford, Fig. (a) shows EF result (average TPR=0.51 
and FPR=1.6E-05) and (b) shows EFM result (average TPR=0.98 and FPR=9.1E-06). 
The averaging is done over 50 runs. Figures (c-d) show the enlargement of Segment 1 in 
(a) and (b), respectively, and Fig. (e-f) show the enlargement of Segment 2 in (a) and (b), 
respectively. Figures (g-h) show the matching, between two small regions in the left and 
right images, of the EF and EFM results, respectively. The average GQ ratios are 1.042 
and 1.0630 for EF estimated models, and 1.0102 and 1.0079 for EFM estimated models. 
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Figure 9: First column shows left images of the examples, second and third columns show 
the EF and EFM results, respectively. The average increase, over 50 runs, in the number 
of matches found by EFM in comparison to EF for the examples shown above (top to 
bottom) is 0.76, 10.53, 3.33, 0.44, 0.6 and 0.68, respectively. 
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GQ 


ROC 


median 


mean 


variance 


TP 


FP 


TPR 


FPR 


EFM 


1.0048 


1.0074 


4.00E-06 


824 


18.08 


0.98 


2.30E-06 


EF 
SBR=0.6 


1.0386 


1.0475 


1.00E-03 


602 


31.66 


0.72 


4.10E-06 


EF 
SBR=0.7 


1.0415 


1.0519 


1.30E-03 


652 


41.14 


0.78 


5.20E-06 


EF 
SBR=0.8 


1.0460 


1.0521 


8.00E-04 


691 


51.60 


0.82 


6.60E-06 



EFM 


1.0183 


1.0194 


1.00E-06 


501 


26.96 


0.97 


3.10E-06 


EF 
SBR=0.6 


1.1742 


1.3031 


1.71E-01 


94 


19.18 


0.18 


2.20E-06 


EF 
SBR=0.7 


1.1989 


1.3012 


8.64E-02 


171 


33.22 


0.33 


3.80E-06 


EF 
SBR=0.8 


1.0806 


1.2594 


1.09E-00 


256 


49.20 


0.49 


5.6E-06 



EFM 


1.0523 


1.0698 


2.00E-03 


300 


15.72 


0.96 


1.70E-06 


EF 
SBR=0.6 


2.6412 


2.6413 


1.30E-06 


9 


2 


0.03 


2.20E-07 


EF 
SBR=0.7 


1.8993 


2.2440 


1.22E-00 


19 


5.48 


0.06 


5.90E-07 


EF 
SBR=0.8 


2.4915 


3.8799 


9.31E-00 


36 


13.04 


0.12 


1.40E-06 



EFM 


1.0102 


1.0102 


1.90E-09 


656 


36.49 


0.98 


9.10E-06 


1.0046 


1.0079 


1.20E-05 


EF 
SBR=0.6 


1.0625 


1.0681 


7.00E-04 


258 


45.48 


0.38 


1.10E-05 


1.0397 


1.0514 


1.80E-03 


EF 
SBR=0.7 


1.0383 


1.0420 


2.00E-04 


344 


63.04 


0.52 


1.60E-05 


1.0427 


1.0630 


2.20E-03 


EF 
SBR=0.8 


1.0218 


1.0243 


3.00E-04 


431 


76.18 


0.64 


1.90E-05 


1.0447 


1.0905 


9.80E-03 



Table 1: In the case of a single model and increasing viewpoint (Graphite VGG Oxford), 
the first three blocks show the average, over 50 runs, ROC attributes and GQ of EFM 
and EF (using different SBR ratios). The EFM and EF results were comparable for small 
viewpoint but as the viewpoint increases EFM model estimates becomes more reliable 
in comparison to EF estimates. The last block shows the results for a multi-model case 
(Metron College VGG Oxford). In both cases, EFM achieved near optimal matching. 



Tech.Rep.#755, CS dept., UWO, London, Canada, ISBN: 978-0-7714-2980-4 



23 



VI 





GQ 


ROC 


median 


mean 


variance 


TP 


FP 


TPR 


FPR 


EFM 


1.0031 


1.0040 


6.3E-6 


529 


34.20 


0.95 


8.5E-06 


1.0225 


1.1251 


0.0207 


EF 
SBR=0.7 


1.1999 


1.2470 


0.0427 


241 


33.12 


0.43 


8.3E-06 


1.2384 


1.2750 


0.0331 



VI 


EFM 


1.0102 


1.0102 


1.9E-9 


656 


36.489 


0.98 


9.1E-06 


1.0046 


1.0079 


1.2E-5 




EF 


1.0383 


1.0427 


0.0002 


344 


63.04 


0.52 


1.6E-05 




SBR=0.7 


1.0427 


1.0630 


0.0022 



co 
VI 


EFM 


1.0084 


1.0083 


0.4E-9 


720 


45 


0.99 


1.1E-05 


1.0046 


1.0079 


1.2E-5 




EF 


1.0372 


1.0347 


2.5E-5 


369 


65.34 


0.51 


1.6E-05 




SBR=0.7 


1.0500 


1.0589 


0.0015 



Table 2: shows the effect of the fitting threshold T (used in computing ground truth, 
EF and EFM) on the average quality of estimated models and ROC attributes, over 50 
runs. The performance for both EF and EFM in the case of T < 2 is better than the 
case of T < 1 because the threshold in the later case is underestimated — see GQ variance. 
Furthermore, the average increase in the model estimate qualityfor EF and EFM between 
T < 3 and T < 2 is not as significant as the average increase between T < 2 and T < 1. 
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A Lemma 13.11 Proof 



The following proof assumes that there exists a feasible solution for GAP with a 
finite objective value. GAP is unfeasible when ^ |J>|, e.g. \J-{\ < \J- r \ and in 
that case adding (|J>| — \ J~i\) dummy features with a fixed matching penalty T to 
Ti will ensure the GAP feasibility. The objective value of a feasible GAP solution 
is guaranteed to be finite when GAP is solved over any set of models and an outlier 
model (f) with a fixed matching penalty T for all possible pairs of matched features. 
We will prove a more general theorem than Lemma |3.1 Lemma 3A_ is a derivative 
of Theorem IA.1I 

Theorem A.l "There exists an optimal solution, with an objective value k* , of a 
GAP instance if and only if there exists a valid MCMF F* over Q* , of the GAP 
instance, with cost(F*) = k* ". 

Proof Assume that there exists a GAP optimal solution M.*j with an objective 
value k* . If there exists a valid MCMF F over Q* with cost(F) = k such that 
k < k* then we can construct a feasible GAP solution M.f where Aif = {x pq h = 



F(n p h,n q h) | p G Fi,q G J- r ,h G £}. Using Corollary A. 2 we can deduce that the 
objective value of the constructed GAP solution Jiif is equal to cost(F). Now we 
prove that the constructed solution Jiif is feasible by showing that A4f can not 
be unfeasible, i.e. one or more of the constraints (JsJ) can not be violated in the 
constructed solution. Constraints (|8|) are violated when 



I a feature p G T\ is not assigned to any feature. 
That means the MCMF F used to construct M. f does not saturate Q* and this 
is a contradiction to our assumption that F is a MCMF. Notice edges (s, n p ) 
for p G Ti must be saturated as in the worst case scenario p will be matched to 
another feature through the outlier model for a fixed cost penalty T . 

II a feature p G T\ is assigned to more than one feature in T r , e.g. q\ and q2- 
If there exist two models h and £ such that x pqi h=l and x pq2 e=l then F(n p h, n qi h)= 
1 and F(n p £,n q2 i) = l must be true by construction of Aif. By construction of 
Q*, n p h and n p £ acquire their flow from n p , and n p could only push out one unit 
of flow. Therefore, for F(n ph , n qih ) = 1 and F(n p £,n q2 i) = 1 to be true n p must 
push two units of flow and that contradicts our assumption that F is a valid 
flow over Q* . 

Ill a feature q G T r is assigned to a zero or more than one feature in J-j. 

We could show that scenario could not happen for Aif by reversing the roles of 
p and q in I and II. 
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IV a matched pair of features p and q are assinged to more than one model, e.g. hi 
and hi. 

If x pqhl =l and 2^2=1 then F(n phl , n qhl ) = 1 and F(n ph2 , n qfl2 ) = 1 must be 
true by construction of M.f. By construction of Q*, n p h x and n p h 2 acquire two 
units flow from n p while n p could only push out one unit of flow. Therefore, for 
F(n phl , n qhl )=l and F(n ph2 , n qh2 ) = l to be true n p must push our two units of 
flow and that contradicts our assumption that F is a valid flow over Q* . 

Finally, if such a solution Aif exist then M*^ is not optimal as k* will be bigger 
than k and that contradicts our main assumption that A^J- an optimal GAP solution, 
i.e. k* is the lowest possible objective value. 

Assume that F* is a valid MCMF over Q* with costiF*) = k*. If there exists an 
feasible solution Aif for which the objective value is k < k* then we can consturct 
a valid MCMF F where 

I G T r where x pqh = 1 

Vp g J 7 /, At g £ 

herwise 

Vp £ Ji,g 6 J>,/i G £ 

> G J} where = 1 

Vg G J>, h E C 

herwise 

Vg G J>. 



F(s,n p ) 




F(n p , n ph ) 




F(n ph ,n qh ) 


Xpqh 


F(n qh ,n q ) 


-{: 


F{n q ,t) 


= 1 



Using Corollary A. 2 we can deduce that the cost(F) = k. 

No we will prove that the constructed flow F is a valid MCMF. A flow is considered 
valid if it satisfies the capacity and conservation of flow constraints over Q*. F satisfies 
the capacity constraints by construction of F — the flow through any edge is either 1 
or while all edge capacities are 1. Furthermore, F was constructed in a way that 
preserves the flow with in Q* . That is, if there is a flow through edge (n p h,n q h) we 
create a flow from s to n p h and from n q h to t along the following paths {s, n p , n p h} and 
{ n qh,n q ,t}, respectively. Therefore, the conservation flow is preserved at n p h,n q h,n p 
and n q . Notice that there can not exist n phl and n pfl2 both creating along {s, n p ,phi} 
and {s, n p , ph 2 } as in this case Ai f will unfesaible which is a contradiction. Moreover, 
the amount of flow going in to Q* through s is \ J-i\ and the amount of flow going out 
of Q* through t is \JF r \ and since Aif is a feasible GAP solution, by definition, then 
\J~i | must be equal to |J>|. Thus, the conservation of flow constraint is preserved at 
s and t, and Q* is saturated. Thus, the constructed flow F is a valid MCMF flow. 

Finally, if such a solution F exist then F* is not MCMF as cost(F) < cost(F*) 
and that contradicts our main assumption that F* is a MCMF over Q* . 
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Corollary A. 2 For a valid F over Q* and a GAP solution Ai f where A4 / = {x pq h = 
F(n p h,n q h) | p G Ti,q G J- r ,h G C}, the objective value of the GAP solution M.f is 
equal to cost(F). 

Proof 

cost(F) = c{v , w) ■ F(t>, w) by definition of flow cost 

= c(n p h, n q h) ■ F(n p h, n q h) by construction of <7*, other edge costs are 

{rtph,n q h)££ 

= D P q(9h) ' F(n p ft, riq/j) by definition of c over £ 

(™ P / l .'Vi)e£ 

= D pq {9h) ■ Xpqh by condition in Corollary A.2 



ge.Fr- 

B Total Modularity Proof 

Our generalization of the assignment problem for solving the matching problem over 
a set of models C such that \C\ > 1 could be formulated as integer linear program 

GAP : argmin E E E D pq (d h )x pqh 

s.t. £ Eaw» = l VgGJ- r (12) 

E E = 1 G J 7 / (13) 

Xpg ft g {0, 1} Vh e £, p e T h q g T r . 

The rest of this section proves that GAP is an integral linear program, that is, its 
LP relaxation is guaranteed to have an integer solution. 

Let us denote the coefficient matrix and the right hand side vector of equa- 



tions (12) and (13) by A and b, respectively. It is known [26] that a linear program 
with constraints Ax = b is integral for any objective function as long as b is integer, 
which is true in our case, and matrix A is totally unimodular. It remains to prove 
that A is totally unimodular. 



Lemma B.l Coefficient matrix A of GAP's linear constraints is totally unimodular. 
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Proof The coefficient matrix A of GAP has a special structure that facilities its 
proof of total unimodularity. Let us assume, without loss of generality, that the 
number of features on the left and right images is n. Then the coefficients matrix in 
case L = {h} could be written as follows 



x llh %lnh 



x nlh x nnh 



A' = 


( 1 • 


. 1 




1 . 


\ 

. 1 




1 




1 




1 






V 


1 




1 




1 / 



eqs (12) 



eqs (13) 



In case \C\ > 1 then coefficients matrix A could be written as follows 



A 



1 

(A' 



2 

A'\ 



A'). 



Heller and Tompkins |27j showed that in order to prove that A is totally unimod- 
ular it is sufficient to prove that the following three conditions are satisfied by the 
coefficient matrix: 

I Every entry of the coefficient matrix is either 0, +1, or -1. 



This condition is satisfied for A by construction, see equations (12) and (13). 

II Every column of the coefficient matrix contains at most two non-zero entries. 
Each column in A corresponds to a unique decision variable, for example x pq h- 



Note that variable x pq h appears only once in linear equations (12) and once in 
linear equations (13). Therefore, variable x pq h appears twice in A. That is, the 
column corresponding to x pq h has exactly two non-zero entries. 

Ill There exists a two set partitioning, say l\ and I2, for the rows of the coefficients 
matrix such that if two non-zero entries in any column have the same sign then 
these two rows are in different sets. And, if the non-zero entries have different 
signs then these two rows belong to the same set. 

Notice that A' satisfies condition III by setting I\ and I2 to the rows of (12) 



and (13), respectively. Also, the coefficients matrix A in case of more than one 



model is simply the horizonal concatenation of the coefficients matrix A to itself 
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|£| times. Thus, the constrains added over the two disjoint sets Ii and I 2 , that 
satisfy condition III over A , by repeating its columns are redundant. Finally, 
condition III will be satisfied by A by the same row partitioning that would 
satisfy condition III for A . 
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